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Lavery splines are examined in the uni- 
variate and bivariate cases. In both 
instances relaxation based algorithms for 
approximate calculation of Lavery splines 
are proposed. Following previous work 
Gilsinn, et al. [7] addressing the bivariate 
case, a rotationally invariant functional is 
assumed. The version of bivariate splines 
proposed in this paper also aims at irregu- 
larly spaced data and uses Hseih-Clough- 
Tocher elements based on the triangulated 
irregular network (TIN) concept. In this 
paper, the univariate case, however, is 
mvestigated in greater detail so as to ftir- 
ther the understanding of the bivariate 
case. 
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1. Introduction 

Over the last decade the task of representing very 
large "point clouds" by meshed surfaces has arisen in 
many applications. Such point clouds may have been 
generated from terrain surveys or from LADAR (LAser 
Direction And Ranging) scans of objects. Work has 
centered around the concept and application of 
Triangulated Irregular Network (TIN). Here the "foot 
prints" iXi,yi) of data points (x^^y^, z^ are "triangulated", 
that is, the area of interest is tiled by non-overlapping 
triangles defining a piecewise triangular "TIN-surface" 
in space as each triangle is lifted according to the ele- 
vations z, at its comers as indicated in Fig. 1, Most 
often, the triangulation used is the "Delaunay" triangu- 
lation, characterized by the "empty circle property" in 
which the circumcircle of each triangle does not con- 
tain any triangle comers in its interior (see Fig. 2), 
Figure 3 features such a Delaunay triangulation, under- 




Fig, 1, TIN meshing: Triangulated surface over a triangulation in the 
footprint plane. 



lying a TIN, created from a LADAR scan of a construc- 
tion site on NIST grounds. The constmction of the sur- 
face from adjacent patches or "elements" is basically a 
"Finite Elemenf technique. 

TIN constructs utilize the actual data points directly 
even if they are distributed irregularly rather than refer- 
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Fig, 2, Delaunay triangulation: Circumcircles of triangles do not 
contain vertices in their interior. 




Fig, 3, A TIN created from the LADAR scan of terrain. 

ring to approximations by regular rectangular grids. 
They adapt naturally to disparate data densities. 
Importantly they support strategic selection of data 
points to be used as the support points for a TIN sur- 
face. 

Conceptually TINs define which points may be con- 
sidered "neighbors" of each other. This proximity con- 
cept based on a direct neighbor criterion differs from 
concepts based on a fixed distance cutoff in that it 
adjusts automatically to differing data densities such as 
in data collected by ground based LADAR. Figure 3 
exemplifies the dramatic differences in densities 
encountered in ground-based LADAR scans. Also, 
because of the neighbor relation provided by a TIN, the 



utility provided by that construct is not restricted to sur- 
face generation, but extends to data editing and analy- 
sis. 

The TIN surface, as defined according to Fig. 1, is 
C^, that is, it is continuous but usually not smooth: it is 
typically not differentiable along triangle edges, since 
the plane of any of the spatial triangles constituting the 
surface generally differs from the planes of its adjacent 
triangles. The TIN concept of triangulation-based sur- 
face generation is frequently understood to imply the 
use of such planar "elements" resulting in a non- 
smooth, piecewise linear TIN surface. It should be 
emphasized, however, that the TIN concept also sup- 
ports the use of non planar, that is, curved elements for 
a smooth or C^ surface. Planar elements offer some 
advantages besides simplicity and a certain robustness 
to be discussed below. Visualization still requires piece- 
wise linear representations in order to delineate hidden 
surfaces. Volume computations similarly tend to be 
either grid-based or based on piecewise linear surfaces 
for ease of computation. In those cases, a smooth sur- 
face would have to be discretized or approximated by a 
piecewise linear surface. This then begs the question, 
why not use planar elements in the first place? 

Without any doubt, however, there are many 
instances, in which a smooth surface representation 
would provide definite advantages: 

• More accuracy could be achieved without increasing 
memory requirements because curved elements 
encode more geometric information than planar ele- 
ments. 

• Clusters of coplanar points are avoided when sam- 
pling the surface at a discrete set of foot prints such 
as regular grid points. 

• When simulating movement over a terrain surface, a 
smooth ride is preferable. A ride over a piecewise lin- 
ear surface is necessarily bumpy at transitions from 
one triangle to another. 

• In many applications, several point clouds which 
were collected with reference to different coordinate 
systems need to be "registered", that is, combined 
within a common coordinate system. Some common 
registration methods, such as point-to-surface itera- 
tive closest point (ICP) proposed by Besl [2], require 
minimizing the deviations of points in one point 
cloud from the surface representation of another 
point cloud. This minimization process should work 
much better if that surface were smooth. 

In computer aided design (CAD), design-defined 
objects are routinely represented by surfaces that are C^ 
or even C^. Terrain representation, however, poses a 
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different challenge: The location of "crease lines" or 
"break lines", along which the actual terrain surfaces 
are not differentiable, are usually not known ahead of 
time, whereas in a CAD environment break lines tend 
to be specified as part of the design. Unspecified break 
lines, on the other hand, along with actual verticalities, 
tend to give rise to spurious oscillation and, what may 
be called, "Gibbs phenomena" in analogy to the phe- 
nomenon well known from the theory of Fourier series. 
TIN surfaces based on planar elements, on the other 
hand, are more robust, and can be constructed so as to 
automatically represent and, if necessary, report break 
lines. Susceptibility to spurious oscillations and Gibbs 
phenomena are one of the reasons why the terrain mod- 
eling community has been slow to accept smooth sur- 
faces. There has thus been a long quest for "non-oscil- 
latory splines" which would obviate this pesky conun- 
drum. 

In 1994, Lavery [11,12,13] (see also Gilsinn and 
Lavery [8,9,10]) proposed successful paradigms for 
univariate as well as bivariate non-oscillatory splines, 
which could be used for representing 2D or 3D data 
sets, respectively. Lavery introduced the term "L^ 
splines" for his brand of nonoscillatory splines. The 
term L^ splines is, however, frequently misinterpreted 
as minimizing an Ly measure-of-fit when approximat- 
ing a point set by, say, classical splines. For this reason, 
we prefer the term "Lavery splines". 

Classical splines are characterized by their minimiz- 
ing energy functionals. Lavery splines, on the other 
hand, minimize different functionals. In the bivariate 
case, especially, the computational effort of minimizing 
these functionals, however, exceeds the effort required 
by the classical approach by an order of magnitude. For 
the bivariate case we have, therefore, considered an 
approximation to the calculation of Lavery splines, 
along with a prior modification of the functional pro- 
posed by Lavery for the bivariate case. This modified 
bivariate functional is still an extension of Lavery 's 
functional for univariate functions, but it extends a dif- 
ferent aspect of the latter. It also is invariant under pla- 
nar rotations of the coordinate system, which Lavery's 
bivariate functional is not. We will present preliminary 
results for both univariate and bivariate nonoscillatory 
splines. 

The main thrust of this paper, however, remains to 
illustrate the utility and performance of bivariate non- 
oscillatory splines, building on the previous study by 
Gilsinn et al. [7] , and also on an early terrain modeling 
study by Mandel et al. [18]. Both studies employ TIN 
techniques in conjunction with the "reduced Hsieh- 
Clough-Tocher (rHCT)" element, the approach was 



pioneered by Lawson [14,15] and is followed in this 
work. It involves specifying elevations z, and two 
slopes z^ , z^y at each triangle comer or vertex v, = {x, , y^ 
in a triangulation, and filling in, at each triangle, the 
thereby defined rHCT elements, results in a smooth 
surface over an entire TIN. 

The classical spline approach to, say, interpolating 
the elevations at vertices of a triangulation would be to 
prescribe the elevations z,, and select the remaining 
parameters, namely the partial slopes z^, z^y, so as to 
minimize some energy functional. This was indeed the 
approach taken in the work by Mandel et al. [18] The 
task to be accomplished there was to represent terrain 
given by digitized contour lines. The elevations at the 
data points were thus determined by the elevation of the 
contour line on which the data point was located. In the 
end, elevations were to be evaluated at the points of a 
square 900 x 900 grid. At the time, — 1985 — visuali- 
zation procedures that are now routine were not yet 
commonly available. A display of the results had to 
wait another year. So assessing the quality of the repre- 
sentation was restricted to manual spot checking. It 
thus took until a day before a scheduled presentation 
that, to our dismay, huge oscillations were detected. 
Fortunately, it turned out that it was not due to a prob- 
lem with the method but was caused by an error in the 
data that had been provided: the elevation of a single 
contour line had been tagged 100 feet too high, causing 
the disruption. In this instance, the oscillations proved 
to be actually beneficial in that they uncovered data 
errors. 

In Sec. 2 we discuss aspects of univariate Lavery 
splines in order to shed light on related issues for 
bivariate splines, which are addressed in Sec. 3. The 
algorithm proposed here for bivariate non-oscillatory 
splines requires solving large sparse systems of linear 
equations. Experience was gathered concerning the 
performance of the well known Gauss-Seidel algo- 
rithm. 



2. Univariate Spline Interpolation 

While the emphasis of this work is on bivariate 
splines for surface generation, an examination of uni- 
variate splines is taken up in order to highlight some of 
the issues pertaining to splines in general. In the uni- 
variate case, cubic "spline functions" are most com- 
monly used and are considered here. They form a linear 
space 
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of piecewise cubic C^ functions f{x) defined locally 
over intervals between "knots" 

Xo<Xi<... <x„, 

that is, they consist of cubic polynomials 

Mx),XG [X,._i,xj, /= 1,..., 77. 

Adjacent cubic polynomials are required to assume the 
same values y^ at common interior knots, 

3^/=y^(^/)=Ai(^/). 

This ensures continuity of the complete spline function 
f(x) over the entire interval 

In addition, the polynomials are to assume the same 
slopes 

The spline functions /(x) are thus continuously differ- 
entiable, that is, they belong to class C\ In what fol- 
lows, the linear spaces 

of first and second derivatives of spline functions are 
also considered, in spite of the fact that, at common 
knots, the second derivatives of adjacent cubic polyno- 
mials may not agree, so that the spline functions/(x) e 
F are generally not twice differentiable at such knots. 
However, they are twice differentiable everywhere but 
on this set of measure zero. For the purposes of integra- 
tion below, it does not matter that the function f'\x) 
may not be defined for those arguments. 

Each of the constituent cubic polynomials f^x) is 
uniquely determined by the values yi_i, y^ at the knots 
Xy_i, Xi and the slopes m,_i, m, at those locations (Fig. 4), 
in fact, the polynomial is linear in these parameters. 
The entire function/(x) is thus uniquely determined by 
its values and slopes at the knots, and it, too, depends 
linearly on these parameters, so that the space F is iso- 
morphic to the 2(«+l)-dimensional vector space of val- 
ues;^, and slopes m^, i = 0, ..., n. 

We now turn to the task of interpolation. Here the 
values y^ at the knots x^ are fixed and specified. Given a 
particular specification of values yi, a corresponding 
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Fig. 4. Hermite cubic polynomial determined by coordinates and 
slopes at end points. 



"interpolating spline function" depends only on the 
parameters mf. 

f(x)=f(mQ,m^, ..., m„;x). 

Collectively, these functions define affine manifolds or 
flats, 

within the linear spaces F, F\ F'\ respectively. That is, 
if 

m.= , z = 0,l,...,«, 

then it follows that correspondingly 

f(x) = f(m^,m,,,„,m^;x) 



The question then becomes, how to select slopes m^ so 
as to achieve a "satisfactory" interpolation. That selec- 
tion is generally made by minimizing a functional on 
the affine space S defined by an integral. In this work, 
we reserve the term "spline" — as opposed to "spline 
function", for the results of such a minimization. 

2.1 Paradigms for Univariate Splines 

"Classical splines" are uniquely defined as those 
interpolating cubic spline functions which are (i) C^, 
that is, twice differentiable, such that 



(1) 



fI'{x^=U{W,i=\,...,n-\ 



holds at interior knots, and for which (ii) the second 
derivative vanishes 



(2) 



f\x,)=r{x„)=o 
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at the two exterior knots. Holladay proved early on (see 
Ahlberg [1]) that classical splines are also defined as 
the unique minimizers of the "thin beam energy" 



(3) 






over the affine space S of all C^ interpolating spline 
functions. 

Condition (2) is familiar to structural engineers as 
the vanishing of the second derivative at the "free end" 
of a beam. It is, however, remarkable that this energy 
minimization enforces also a higher level of compati- 
bility across knots (I) so that the minimizing C^ spline 
functions are, in fact, in class C?. On the other hand, 
this stiffness contributes to the tendency of classical 
splines to produce spurious oscillations, Gibbs phe- 
nomena and undesirable inflections. 

Several attempts have been made to avoid these 
problems. Taking a clue from mechanics, Schweikert 
[20] and Cline [3] have introduced "tension splines", 
where the arclength of the spline function is made part 
of the defining minimization. Reinsch (see Stoer and 
Bulirsch [21]) moved to "exponendal splines" by 
adding a multiple of the square of first derivatives to 
the integrand in (3). Those efforts were only partially 
successful. Drawbacks include the need to specify an 
additional parameter in order to balance conflicting 
minimization requirements, and the fact that these tech- 
niques are not readily generalized to the bivariate case. 

Lavery splines, on the other hand, appear to avoid 
such shortcomings. In Figs. 5 and 6 we compare Lavery 
splines against classic splines. We are particularly 
impressed with the performance of Lavery splines for 
the example in Fig. 8, as compared to the example in 
Fig. 7. What is the secret of Lavery splines? How are 
they defined as opposed to classical splines? 

Lavery defines his nonoscillatory splines, in essence, 
as minimizing the integral of the absolute value rather 
than the square of the second derivative of a spline 
function: 



(4) 
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Fig. 5. Lavery spline vs. classic spline: Example 1. 
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over the given affine space S of interpolating spline 
functions. 
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Fig. 7. Classic splines produce unnecessary variations in curvature. 



2.2 Expressing Holladay and Lavery Integrals 

For the purposes of this paper, we will refer to the 
integrals (3) and (4) as the "Holladay integral" and the 
"Lavery integral", respectively. The goal of this section 
is to derive expressions for the values of these integrals 
in terms of the slopes rui, i = 1, ..., n, at the knots of the 
spline function /(x) under consideration. Both are the 
sums of the corresponding integrals of the second 
derivatives of the individual cubic polynomials yj(x), 
which constitute the spline function/(x). 

As pointed out in Sec. 2.1, each such cubic polyno- 
mial is uniquely defined by its end points (X/_i, j/_i), (x^, 
y^ and its end slopes m,_i, m,. The various formulas 
describing this polynomial are commonly referred to as 
"Hermite" formulas. For the versions used here, we 
introduce the quantities 



and 



M, 



A,- = Xr 



yi-yt- 



yi-yt- 




-6 -4 -2 a 2 4 

Fig. 8. Lavery spline has smooth curvature transitions. 



where M^ represents the slope of the straight line 
between end points. Instead of referring to the variable 
X directly, the following formulas are in terms of the 
weights 



(5) 



A^.=A^.(x) = -^ 



ix.=ji.{x) = 



^i-^i-l A- 



where A, + fii= 1, and A,, jti, > 0, for x in the interval 
\pCi-u -^J - Such weights are often referred to as 
"barycentric coordinates" or, in the bivariate case, as 
"triangle coordinates". With these conventions, we 
find, for instance, 

(6) fix) = 

V/-1 + liyi + ^]liiijni-i - M-)A,- - X,ji]{m, - M-)A, . 

Furthermore, by definition (5), 

dr = -A,dA, = +A4jU,., 

the chain rule yields, using (A, + ji^j^ = 1, 
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fix) = 6X,^Mi + (A/ - 2A,A>/-i + (// • - 2 V,)m, 

An alternate expression for the first derivative is readi- 
ly derived: 



(7) fix) = A,m^i + //,m, + A^A ■ 
Here the quantity 

(8) D. = 6M. -3m._i -3m. =6\ M. '^ '- 



vanishes if and only if the polynomial /(x) has degree 
less than three. Because by (5) 



f^ 



Y^ 



X.fi. = 



-Di is seen as the lead coefficient of/'(x) as expressed 
inx. Consequently -1/3A is the lead coefficient offjix), 
Di < indicates that the function is concave up to its 
inflection point and convex thereafter. Conversely D^ > 
indicates that convexity precedes concavity in the 
direction of the x-axis. The quantity D, will play a 
major role in what follows. The same is true for the 
quantities t/,, V^ in the expression 



(9) 
where 

(10) 



y; (x)=— a,i7,-/i,j;), 






(11) Eif,) = 

— [m._j +m._^m. +m. -3M.(m.^ +m.)+3M. ]. 
A, 

The full Holladay integral is the inhomogeneous quad- 
ratic function of the variables m^, arrived at by adding 
the energies E(f^ of all partial functions fix). 
Introducing the vector of slopes 

m=(mQ,m^, ...,my, 

we have in matrix notation 

E(f) = m'Hm - MM^m + 12c, 

where 



H = 



_4_ Jl_ 

A, A, 

J^ J. y_Jl_ Jl_ 

Ai Ai A2 A2 

A2 A2 A3 A3 

2 2 I 2 J 

A„_i A^i A„ A, 

2 _4 

A.. A, 



M = 



M, 






A, 






A2 


+ 


Ml 


A3 


+ 
+ 


M, 
A2 


M„ 


+ 


M„_ 
~- 

"a^ 



The coefficients t/,, V^ thus relate to the two-sided sec- 
ond derivatives of ^(x) at knots x, . Note that 

and also that, in view of (5), 



:-^, £ ^Mdx = ^, £ fifdx = ^ 
3 ''^M 6 •'^^1 3 



An expression for the thin beam energy of the individ- 
ual polynomial ^(x) is now readily derived: 



and 



.^1 



M 



M„ 



c={-^yH-^y+,-, H-^r- 



Minimizing this expression for the Holladay energy 
integral is to solve the linear system of equations 



(12) 



or 



2Hm=M 
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2 1 3 ,^ 

— nif. H nt = — M, 



{ 



^0 + 



2 2 



A 



+ 

y \ L J 



1 



m + — nu - — M + — My 

K K K 

'2 ^1 



A, " A ' A2 



( 



-m, + 



2 2 



A 



y 2 3 y 



1 3 3 

A3 ^ A, ^ A3 ^ 



1 2 3 ^^ 

A " A A 



The first equation, in fact, may be restated by (10) as 

that is, as the requirement (2) that second derivatives 
vanish at the first knot. The last equation reflects the 
corresponding requirement at the last knot. The second 
equation is equivalent to 

— ^1+— t/, =-/'U)--//(;q)=0, 

A^ ' A2 2 2 

implying second order differentiability at the interior 
knot Xy, The remaining equations similarly enforce 
compatibility of the one-sided second derivatives at the 
remaining interior knots. This confirms the result of 
Holladay 

The reader should note that the system (12) for clas- 
sical splines differs from the linear system mostly 
offered in the current literature. There the second order 
differentiability of the classical splines is already 
assumed and the spline functions are formulated in 
terms of those second order derivatives, n^ =f'\x^, i = 
0, ..., n. However for weighted classical splines, to be 
encountered in Sec. 2.4.1, second order differentiabili- 
ty no longer holds, and a weighted version of the linear 
system (12) needs to be considered. 

2.2.1 Expressing the Lavery Integral 

As to the integral of the absolute value of the second 
derivative, it is readily available in the case that the 
derivative does not change sign inside the subinterval: 



e(y;)=riyr(x)idx=r/(x)dx 

=1 fi(^t)-fi(^t-i) \ = \^ -^i-i I ■ 



This includes the case in which the polynomial is of 
lesser degree than cubic, and thus has constant second 



derivatives. This case is signaled by the vanishing of 
the quantity D^ introduced earlier in (8). 

The function /''(x), however, is a linear function in x 
and, unless constant, changes sign at some location x, 
which also marks the location of the inflection point of 
fi{x). Suppose this location falls into the interior of the 
subinterval: 

X,-_i < X < X;, 

Then the integral has to be calculated as the sum of two 
integrals of linear functions: 



eifd- 



f /(x)dx|+|r/(x)dx. 



The integrals between the absolute value bars are of 
opposite signs, so that the total integral can be written 
as the absolute value of a difference of integrals 



e{fd-- 



J* /(x)dx-j:'/(x)d^. 



These integrals can be separately evaluated as differ- 
ences of slopes. Let 



denote the inflection slope, then 
interior inflection — 



in the case of an 



(13) 



e(/3 = |2m-m^_i-m,-|. 



The quantities Ui, V^ depend linearly on the two slopes 
mi_i. Ml. Conversely, those slopes can be expressed in 
terms oft/,, V^. 

(14) m,_, = M, -hj^ +lf^, m, = M,. +^-U,, -| V, 

The next step is to express the inflection slope in terms 
of Ui, V^. To this end, we express the inflection argu- 
ment by its barycentric weights: 

(15) X = ;[,x,_i + //>„ 

From the expression (10) forffXx) it follows that 






where the denominator D^ = Ui+ V^ = 6M^ - 3m^i -3m/ 
has already been encountered in (8) as a quantity that 
vanishes if and only if the polynomial in question is 
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parabolic or linear. Thus Z), = leads back to the previ- 
ous case of no sign changes by the second derivative. 

The weights Xi, fl^ may now be inserted into the 
expression (7) for the derivative /-(x). This gives 

A A A 

Substituting for m^_i, m, according to (14) yields 

m. =M. +— ^ 



3D, 



as well as 
(16) 



e(f^-- 



t2 , Tz2 ,'i2 , 2^ 



Uf+V, 



M, 



n 



3|AI 3 

in view of (19). In terms of the slopes m^i, m{. 



eifd-- 



1 8M/ - 1 8M. ( m._i + m. ) + 5 mf_i + 8 m._^ m. +5 ?7^^ 



6M, -3 m, 1 -3 m, 



The above expressions for e(/J) are also valid if the 
inflections occur at the ends X/_i, x, of the interval of 
definition, reducing to e{fi) = |m, - m^_i|, in accordance 
with earlier results. 

We are now ready to examine the full Lavery inte- 
gral. At first blush, all that remains to be done is to sum 
over the partial integrals e(/J) in their various forms. We 
will show, however, that many terms of the expressions 
(13) cancel each other out as these partial integrals are 
added together. To this end, we distinguish five sepa- 
rate kinds of polynomials f^x) depending on their 
behavior in the interior of the interval between its 
knots, x,_i < X < X,: 

• "Linear"; herey^"(^) = ^ throughout 
and e(/D = 

• "Convex"; herey^"(^) ^ in the interior of the inter- 
val [x,_i, xj 

and e(fi) = m,- - m^_i 

• "Concave"; herey^"(^) < in the interior of the inter- 
val [x,_i, xj 

and eifi) = m,_i - m,- 

• "Convex-concave"; hcrcf-Xx) > for x < x, f-Xx) < 
for X > X 

and e(/J) = +m- m/_i - m„ also Z), ^ 

• "Concave-convex"; hcrcf-Xx) < for x < x, f-Xx) > 
for X > X 

and e{f^ = -m + m^_i + m,, also D^ < 

The last two categories are the ones with an interior 
inflection point x and inflection slope m . 



The interior inflection points considered so far may 
not be the only inflection points of the spline function 
/(x). Inflections occur also at knots x^, < z < n if a con- 
cave or convex-concave polynomial is followed by a 
convex or convex-concave polynomial and, analogous- 
ly, if a convex or concave-convex polynomial is fol- 
lowed by a concave or concave-convex polynomial. We 
refer to such knots as "inflection knots". To make mat- 
ters more complicated, however, inflection may occur 
along an entire stretch of consecutive linear polynomi- 
als of equal slope, the inflection slope in this case, pro- 
vided there are adjacent nonlinear polynomials at both 
ends of such a stretch exhibiting the same 
convexity/concavity pattern that would cause an inflec- 
tion at a knot. In this case, we choose an arbitrary knot, 
say, the leftmost one in the linear stretch, as the inflec- 
tion knot representing the inflection. 

Clearly, slopes at interior knots that are not inflec- 
tions cancel out as the expressions (13) for the Lavery 
integrals e(/J) for the partial spline functions f^x) are 
added up. See Fig. 9 for an example. This leads to 

Observation A: The Lavery integral of a cubic 
spline function is the absolute value of an alternating 
sum of the inflection slopes and of the end slopes. Let 



mi, m2, ^3, 



mr 



be the sequence of all inflection points identified above, 
sorted from left to right by their location. Then 




]\ri^} 



dx - 



x^ x^ x^ x^ 



- \~ ^^0 - 2 W2 1 -2/^3 +2/^5 - m 5 1 



Fig, 9. Lavery integral as expressed by inflection slopes and end 
slopes. 
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e{f) = 



,+2£(-l)'«,-(-l)^ 



(Note that the indices / of m/ do not refer to the interval 
in which they are located). 

2.3 Properties of Lavery Splines 

In this section, we gather some information about 
Lavery Splines, namely, interpolating spline functions 
in the affine space S which minimize their respective 
Lavery integrals. A first general observation concerns 
the convexity of the Lavery integral. 

The Holladay and Lavery integrals (3) and (4) of a 
piecewise cubic spline function, 

/(x) =/(mo, mi, ..., m„;x) 

are functions of the slope specifications mf, 

E{f) = E(my m2, ..., mj, e(f) = e(my ^2, ..., mj. 

The quadratic function E(f) representing the Holladay 
integral can be shown to be positive definite and, there- 
fore, strictly convex. Its minimum is unique on the 
affine manifold S interpolating spline functions. The 
restriction to S is, of course, necessary as the value of 
E(f) would not change if the spline function/(x) were 
modified by adding a linear function in x. Adding a 
non-zero linear function to the function/(x) would not 
preserve its interpolation property, 

2.3.1 Convexity 

Next, we establish the convexity of the Lavery inte- 
gral. It may be viewed as the extension of the Ly vector 
norm to a norm on the linear space f of second deriv- 
atives of spline functions: 

eif) = iiriii. 

The following generic seminorm properties are easily 
verified for the piecewise linear functions in F\ 

ll/"lli = 0<=>/" = 

l«ll/"llll = l«lll/"lll 

along with the triangle inequality. 



Suppose the two spline functions /^^\/^^ are actual- 
ly two interpolating spline functions and, therefore, in 
the affine space S, Then their mean is again in S and, 
from the triangle inequality, 

„r'+/^^'„.iir' 11+11/^^' II 



In terms of Lavery integrals. 



nf^^ + mf^ 



m^^+nf? 



I 2 ^ J 



which establishes convexity. This leads immediately to 

Observation B: Any positive linear combination 
of Lavery splines for the same interpolation problem, 
— in particular, their mean — , is again a Lavery spline 
for this problem, 

2.3.2 Uniqueness of Inflections 

Contrary to the Holladay functional, which is strict- 
ly convex, the Lavery integral is not. As a result, 
uniqueness does not follow and, in fact, does not hold, 
as an example in Sec. 2,3.4 will show. Such minima of 
a convex function, however, must form a convex set. 

Note that, in general, 

J 2 J 2 ' 

If both/^^^ and/^^ are Lavery splines for the same inter- 
polation problem, then so is their mean, and all three 
functions 



(1)* , Mir 



r^ +/ 



M\r Mif 



return the same optimal value for their Lavery inte- 
grals. Thus equality holds in the above relation. This 
implies that both /^"W ^^^ f^'X^) have the same sign 
pattern: 

P'^\x) > if and only iff^^'Xx) > 0, 

This can be rephrased as 
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Observation C: Two Lavery splines for the same 
optimization problem share essentially the same inflec- 
tions: if one of them has an inflection point at x = x, 
then so has the other unless it is linear at this point. 

2.3.3 Free Ends of Lavery Splines 



Xi(6M^ - 3mQ - 3mi) = 3M^ - nt^- Im^, 

which — for the particular value (17) of >^i — yields 
the corresponding end-slope 



(18) 



m^ 



io-Vio_ 5-A^ 



-M 



-n\. 



Here we will examine the free ends of Lavery 
splines, in particular, the cubic polynomial /i(x) and the 
corresponding first summand e(/i) of the Lavery inte- 
gral e(f) . As the end slope mo may vary freely, it must 
optimize eif) while keeping the slope m^ fixed. This 
fact determines the behavior of Lavery splines at free 
ends. As seen in the previous section, the first summand 
e(f]) is a convex fiinction in the variables mo, m^. If the 
slope m^ is held fixed, e(f)(mQ) is convex as a function 
in mQ alone. If the fixed slope equals the straight-line 
slope, then mQ = mi=M^ obviously represents the opti- 
mal value for mo, since /i(x) in that case is a straight line 
with e(f]) = 0. We suppose, therefore, that m^ ^ My, and 
we examine the case, i\mif{x) has an inflection x in the 
interval between the first two knots, Xq<x<Xi, 

Consider the function 






where m as well as t/j, V^ also depend on the variable 
mo. Clearly, the absolute value of ^(/i) is given by e(f^. 
Note that 






in view of (10), and 



= A,'-A'-6A,A 



^ = -2,-^ = -!^ = -3. 

dm^ ' dm^ dm^ 

Solving the quadratic equation for X, and taking into 
account that < >^ < 1, yields 



(17) 



A, =^±^ = 0.860378, 



/^i=- 



6 
4-VlO 



= 0.139622. 



The value of mo for which this value for X^ is realized 
can be inferred from the general definition (15) of an 
inflection, giving the barycentric coordinate X^ in terms 
of the slopes mo, m^, as follows: 



This value for mo represents a locally unique stationary 
value of ^(/i)(mo). 

Now eifXmo) = would imply e(/i)(mo) = 0, and 
consequently linearity, that is, mo = mi = M^, which has 
been ruled out. By continuity, e(f^{m^ is either always 
positive or always negative, — in other words, either 



or 



eifi){mQ) = -eif^Xmo). 

This implies that the value (18) for mo is also a locally 
unique stationary value of e(/i)(mo). In view of the con- 
vexity of this function, it is also its minimizer. At the 
last end we differentiate 



(19) 






3Z)„ 



and find 

de(f„) _ U'„-V„'-6U„V„ 

Symmetrically, we thus have 



= 6/i'-4/i„-l. 



fin - ^\^ K - fii- 

The equation 

fini^K - 3m,_i - 3m,) = 3M, - Im^y - m, 
thus yields a symmetric relationship to (18): 



M„ 



-m„_ 



(20) m„=- 

This establishes 



Observation D: The free ends of Lavery splines 
are either linear functions or they contain an inflection 
in the interior of their interval of definition. The loca- 
tions Xi, x„ of these inflections are universally given, 
respectively by 
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2+Vio 4-VlO 



A Non-U nique Lnrwery Sp line 



-^. 



4->/l0 2+VlO 



-x^. 



Observation D enables us to determine universal val- 
ues for partial Lavery integrals at the end-intervals of 
Lavery splines. By (15), (16), (18), and again by (20), 
which implies 

10-VlO 



A +^ =\+fin =- 



we find 



_^ IQ-VlQ , ^ I ^^, IO-a/iO i r. 



27 



27 



Substituting for jtiq and m„, respectively in D^ = 6M^ 

Sthq - 3mi and Z)„ = 6M^ - 3m„_i - 3m„, 






giving rise to 

Observation E: A necessary, but far from suffi- 
cient, condition for a spline function f to be a Lavery 
spline is that 

e{f,) = ^{M-\)\M,-n,,\, 

e(/„) = |(>/lO-l)|M„-m„.J. 

2.3.4 Examples of Non-unique Lavery Splines 

In this section, we present an example in which the 
Lavery splines are not unique. Consider the three points 
(Figs. 10, 11, 12): 



(21) 



A = fe3^2) = (+1,-1). 



The associated interpolating cubic spline functions 
f{x) are then defined by their slopes 




Fig. 10. One of three Lavery splines for the same interpolation prob- 
lem. 
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Fig. 11. Second of three Lavery splines for the same interpolation 
problem. 
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Fig. 12. Third of three Lavery splines for the same interpolation 
problem. 



at these points. There are two subintervals with cubic 
polynomials 

flAfiiA 



each of them a fi*ee end. This determines the coordi- 
nates ii, X2 at inflections: 
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. 2 + VlQ . _ 2+VIq 

x^ — - ,x^—-\ - . 



Note that both partial Lavery integrals are end-inte- 
grals. For an interpolating spline function /(x) to be a 
Lavery spline it will be necessary by Observation E that 



For mi = +1, the resulting Lavery spline is the symmet- 
ric image of the previous one: 



/^(x) = +x 



2V1O , 10+2^^ 



x^ +x. 



</;)=-(Vio-i)ii 



■my 



_^ Vio-i , . 

</2)= 1-1-^ 



SO that 



</)=-(VIo-i)(|i-«,|+|-i-«^l). 

Clearly 

(22) -1 <mi< 1 

implies 

\\-m,\ + \-l-m,\ = \i\-m,) + i-\-m,)\=2\m,\<2, 
whereas either mi < -1 or mi > 1 would imply 

|l-mi| + |-l-mi| = 2>2|mi|. 

Thus condition (22) characterizes all Lavery splines for 
the example. 

Consider now any slope mi from -1 through +1. For 
mi = -1, we have by (18), (20), and in view of Mi = 1, 
M2 = -l: 

+I5-2VI0 
m^= ,m2 =-1. 

Symmetrically, we find for mi = +1 that 

-I5+2V1O 

m^ = +1, m^ = . 

These slopes, respectively determine the two extreme 
Lavery splines, because each partial function is at a free 
end, and is optimized according to Observation D in the 
previous section. The two resulting Lavery splines are 
extreme in that each has a straight line segment as a 
partial function. 

Using Hermite's formula (6) and substituting for x, 
we fmd for the choice mi = -1, 

^, , 2V1O 3 10 + 2^/l0 

f^(x) = -x. 



-x^ -X 



Both splines are shown in Figs. 10, 11. The self-sym- 
metric spline from the choice mi = is shown in Fig. 
12. In the latter case, 

io-a/Io 

m^= = -m2, 

and 

.^ ^ Vio 3 5+Vio 2 

'5 5 

This Lavery spline is the mean of the two splines with 
linear free ends. This is an instance of observation B 
about positive linear combinations of Lavery splines in 
Sec. 2.3.1. 



2.4 Computing Univariate Lavery Splines 

We now turn our attention to the computation of 
Lavery splines. The commonly used approach (see 
Lavery [11,12]) is to minimize the Lavery integral in 
discretized form, say, 

(23) j^" I r(x) I dx = x^ii rc^.-i ^^^) i 

where the integrand is sampled in each subinterval [x^_i, 
xj at ki equidistant points. Now 

X. , + A. ^— ^—^ — ^ ^—^ 



= \,x._i+//.,x.\, +//., =1. 



Thus 



l'y'(x)<ix=J^^±\X,,U,+H, 



a^l 



= I-rIl^a(3^. -2'",-. -/«,)+«, (3M -'«,..> -2"?)l 

The Lavery integral is discretized as a sum of absolute 
values of the variables m, of the minimization. 
Minimizing such an expression is a Linear 
Programming problem. Many satisfactory methods for 
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solving it are available, such as the Simplex Method or 
Interior Point Methods. The accuracy of the discretiza- 
tion increases with the number of sample points, but 
computational effort increases accordingly 

For that reason, and also to motivate an analogous 
approach in the bivariate case, we are proposing to min- 
imize a different approximation to the Lavery integral, 
one that takes advantage of the ease of computation 
offered by energy minimization. 

By the mean value theorem of integral calculus, 
there exist arguments u^ such that 

£/(x)^dx = A,.(/(M,)^ / = !,...,«. 

we then propose to approximate the Lavery integral by 
the following Riemann sum: 

(24) fir(x)|dx^£A,.|/(«,.)|. 

In contrast to the approximation (23) by discretization, 
this approximation does not offer the option of further 
refinement, unless the interpolation problem itself is 
changed by adding additional knots and ordinates. In a 
sense, it approximates the Lavery paradigm itself rather 
than the Lavery integral. We still use, however, the term 
"approximate Lavery integral" for our proposed 
approximation (24). 

2.4.1 An Iterative Algorithm for Approximate 
Lavery Splines 

For the purpose of computation, we rewrite the non- 
zero terms in (46) 



At each step / = 0, 1, 2, ..., weights 






^/Il, f:i4- 



dx 



r y;'(x)^ck. 



This expression suggests an iterative approach. 
Starting with the classic spline / ^^\x), a sequence of 
interpolating spline functions is generated in hopes to 
converge towards the approximate Lavery integral (24) 



rxxirxx), ...,r\x), ... 

with associated partial Holladay integrals 



E^'-nn^) 



dx, /=],..., n. 



w] 



,(0 



Eun 



,/ = ],..., 77, 



are introduced. Given the function /^^, the subsequent 
function / ^^^^^ is determined as the solution to the fol- 
lowing minimization problem: 

(25) min£w<"£fy;<«''(^|dx. 

This approach raises the question what to do if E(fl^) 
vanishes as the corresponding weight is then not 
defined? 

Simply ignoring such terms may prematurely lock in 
straight line segments between knots. What first comes 
to mind is to specify a limit £ > and boost lower val- 
ues of E(fP) to this level. A more diligent procedure 
might be to start with all weights at value 1, — the 
weight setting that yields the initial classical spline 
according to Holliday's observation — , and then pro- 
gressively increase the use of weights. Such strategies 
remain to be explored. 

In general, adding up partial Holliday integrals, each 
with weight, say, w^ leads again to an expression of the 
energy of a physical structure: a collection of thin 
beams of different thicknesses given by w^, respective- 
ly, and welded together at knots. Minimizing this ener- 
gy expression requires an adjustment to the linear sys- 
tem (12). 

The weighted energy of the partial spline functions 
^(x) is just the product of the straight Holliday integral 
and the respective weight: 

^^ = >v^(/J). 

The total weighted energy is thus given by 

i i 

This means that in the expression (11) for E{f^, the factor 

^ 

is affected as it is multiplied by Wi. In other words, 1 the 
substitution, 



1 w. 



transforms the linear system (12) for minimizing E{f) 
into the linear system for minimizing E^(f): 



70 



Volume 111, Number 2, March- April 2006 

Journal of Research of the National Institute of Standards and Technology 



(26) 



2iVi 



— Lm zz ^-M, 



r}\ +-^m^ = — -M^ + — -M^ 
A, A, A, 

^2 +^7713 = — -M^ + — ^71/3 
A3 A^ A3 



2w^ 2w. 

2_+ — I 



J 



TV 2w 3w 

A„ A„ A„ 



Note that the weighted classic sphnes are, in general, 
not twice differentiable at the knots. However, due to 
the fact that the first and the last equations have com- 
mon factors Wo? ^m respectively, they are equivalent to 
the corresponding equations in the Holladay system 
(12): free ends thus have zero second derivatives in the 
weighted case, too. Again, this is to be expected from 
Physics. 

Note also that some commonly used methods for 
determining classic splines such as B-splines or linear 
systems formulated in terms of second derivatives at 
knots do not carry over to the weighted case. However, 
the above linear system is still "banded", and many 
excellent methods are known for its solution. To solve 
this system we used the venerable Gauss-Seidel 
method, not just for ease of programming, but also 
because it seems to work for bivariate weighted splines. 
Its advantage lies in the fact that the matrix of the lin- 
ear system need not be changed and can be read, so to 
speak, in sequence. This is important for the very large 
systems likely to arise in the bivariate case. In the uni- 
variate case, the convergence behavior is well under- 
stood (See Varga [22]). Using an iterative method for 
solving the class of linear systems above will result in 
a two-tiered iteration procedure: an "outer" iteration, 
developing new sets of weights, and an "inner" itera- 
tion, solving the resulting linear system. Such proce- 
dures can be "balanced", that is, the inner iteration may 
be terminated at a lower accuracy level during the early 
stages of the outer iteration and may be carried to a 
higher level of accuracy as the outer method approach- 
es convergence. This is an added advantage of an itera- 
tive method for solving the linear systems at hand. 

Note finally, that the approximate Lavery method 
proposed here will definitely not converge to the opti- 
mal Lavery integral. This is because the approximate 
solutions are based on weighted classic splines, and 
therefore have vanishing second derivatives at the free 



ends. The second derivatives of Lavery splines, on the 
other hand, assume their inflections in interiors of the 
end interval (See Observation D). 

Observation F: Unless a free end function of an 
approximate Lavery spline is linear, it disagrees with 
the corresponding end function of the true Lavery 
spline in that the latter has an inflection in the interior 
of at least one end interval, whereas the former 
assumes its inflection at the end knot. 

The justification for introducing the approximate 
Lavery splines concept lies in its computational ease 
and the fact that it appears to retain the anti-oscillatory 
dynamic of the original Lavery concept. In fact, the 
examples in Figs. 5, 6, and 8 were calculated using the 
approximate algorithm outlined in this section. 



3, Bivariate Spline Interpolation 

The development of bivariate splines parallels that of 
univariate ones to a large extent. For both, the first step 
is to specify a linear space F of bivariate spline func- 
tions /(x,j), defined on a polygonally partitioned 
region 



Q 



In view of our emphasis on the TIN framework of tri- 
angulations, however, we choose for our discussion a 
rectangular region Q and an irregular set of "knots" in 

V/ = (x,,j/)e Q, /= 1, ...,77. 

An interpolation problem results, if "elevations" z, are 
specified at the vertices v, = (x,, y^. 

The region Q is then partitioned into triangles 



4. 



The vertices 

of any such triangle 4 are among the specified knots, 
and those are the only knots in this triangle. In this fash- 
ion, a "triangulation" of the region Q is obtained. 
Spline functions are then defined by installing partial 
functions 
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Fig. 13. Reduced Hsieh-Clough-Tocher (rHCT) element, zi^, ziy 
denote the partial slopes zx^, zy^ (27). 



in their respective triangles 4 of the triangulation. 

Following C. Lawson [14,15], we install above each 
triangle 4 a reduced Hsieh-Clough-Tocher (rHTC) ele- 
ment. The reduced Hsieh-Clough-Tocher (rHCT) ele- 
ment, displayed in Fig. 14, provides a means for repre- 
senting smooth TIN surfaces (see Lawson [14,15]). The 
rHCT element is defined as a function over a triangle in 
the X, j-plane, and is fully determined by the elevations 
and the partial slopes or derivatives at triangle comers. 
The triangle is divided into three subtriangles with their 
common vertex at the centroid of the three comers of 
the full triangle. In each of the three subtriangles the 
rHCT function is a bivariate cubic function. Together, 
these functions form a function on the full triangle 
(Fig. 13). This smooth rHCT function is furthermore 
constrained by the linear perpendicularity condition on 
the derivatives orthogonal to the outside boundary 
edges of the triangular element: along each such edge, 
these derivatives are required to be linear functions. It 
follows that the orthogonal derivatives along an edge 
are fully determined by their values at the ends of the 
edges. These values are inherent to both triangles adja- 
cent to the edge so that the orthogonal derivatives coin- 
cide when calculated within each of the two triangles 
independently. The above condition thus provides the 
key property ensuring smoothness of a rHTC surface. 



As will be seen in Sec. 3.1.2, the restriction inherent in 
the linear perpendicularity condition may be more 
severe than commonly expected. 
To this end, "elevations" 



and "partial slopes" 

(27) zXi,zyi,i= 1, ..., /2, 

are prescribed at the knots v^. As the vertices of each tri- 
angle 4 are among the specified knots, the correspon- 
ding elevations 



\l'\2'\3' 



and partial slopes 



ZX; , ZV: , ZX; , ZV : , ZX- , ZV : . 



furnish the parameters necessary for defining the rHTC 
elementy^(x, j). As pointed out above in the current sec- 
tion, each specification of elevations and partial slopes 
at the knots yields a C^ function/(x, j), which also rep- 
resents a smooth surface over the region of definition 
Q. All such functions, — relating to the given triangu- 
lation of Q — , constitute a linear space of "spline func- 
tions" and their respective derivative linear spaces 

p p p p p p 

If a specific set of elevations z^ is to be interpolated, 
then those spline functions which meet these elevations 
constitute an affine space of "interpolating spline func- 
tions" and its derivative affine spaces 

^ ^ ^ ^ ^ ^ 

^5 ^x> ^y> ^xxi ^xyi ^yy' 

Note, again, that second derivatives are not fully 
defined for rHTC elements. Indeed, along edges sepa- 
rating two subtriangles of an rHTC element as well as 
along edges separating two triangles of the triangula- 
tion there may be two different values in contention 
depending from which side a point on such an edge is 
approached. For purposes of integration, this is not an 
issue since the discrepancies are restricted to a set of 
measure zero, so that the integrals considered below are 
still well defined. 
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3.1 Paradigms for Non-Oscillatory Bivariate 
Splines 

For the purpose of interpolation, the elevations z, at 
knots are specified, and the task at hand is to determine 
the partial slopes zx^, zy^ so as to arrive at a "best" inter- 
polating spline function. Such an optimal spline func- 
tion is then referred to as a "(bivariate) spline". 
Obviously, there are various kinds of splines depending 
on the choice of the linear space of spline functions and 
the choice of the optimization criterion. 

Paralleling the classical univariate approach, "thin 
plate" minimization exemplifies the "classical" 
approach, and the resulting spline will be referred to as 
the "classical spline", the functional requiring mini- 
mization being 

(28) £(/) = Jj/>2/J+/;J,)dxdj.. 

An approach to computing the surface energy formu- 
la has been announced by McClain and Witzgall [16] 
for rHCT element. A newer version of the report is cur- 
rently being prepared by McClain, Witzgall, and 
Gilsinn [17]. 

The authors believe that, contrary to univariate clas- 
sical splines, bivariate classical splines as defined 
above are, in general, not in class C^, that is, twice dif- 
ferentiable. Nevertheless they are "stiff' and thus sub- 
ject to the dreaded spurious oscillations and Gibbs phe- 
nomena. This can be seen along edges of buildings in 
the urban data set of Baltimore, MD, given in Fig. 14. 
Lavery and Gilsinn [13] were able to address this prob- 
lem by extending his paradigm for univariate nonoscil- 



BaltimoreBOO x GOO Cells 




latory splines to the bivariate case by minimizing the 
functional 

(29) eif) = la /^ I +2 1 4 I +1 4 \)dxdy. 

Again, the squares of second derivatives in an energy 
integral, — here E(f) as defined in (28) — , are replaced 
by absolute values. We will refer to the so defined 
splines also a "(bivariate) Lavery splines". 

In view of the non-uniqueness result for univariate 
Lavery splines, the bivariate Lavery splines are expect- 
ed not to be unique either, as will be discussed below. 
Therefore, the Lavery paradigm typically includes a 
regulatory term such as adding a small multiple of 
/ \x, yY. The resulting bivariate Lavery splines have 
produced oscillation-free interpolations and approxi- 
mations in a large array of applications. They are com- 
puted by discretizations leading to very large linear 
programming problems. 

For large surface generation tasks based on interpo- 
lation, the resulting computational effort may be pro- 
hibitive. Moreover, the bivariate Lavery integral (29) 
functional is not invariant under rotation of the coordi- 
nate system. The authors have therefore proposed an 
alternate extension of the univariate Lavery paradigm 
to the bivariate case: 

(30) e{f) = Ufi+2f^+fl)<ixdy. 

Note, that this paradigm is indeed an extension of the 
univariate Lavery paradigm. Indeed, 



jir(x)idx=jV7Wdx, 



Fig. 14. Urban scene (Baltimore): Gibbs phenomena are visible. 



SO that the integrand of the univariate Lavery integral 
may also be interpreted as the squareroot of an energy, 
— in this case, the elastic energy of a thin beam. 
Replacing that energy expression by that for a thin 
plate, yields the alternate extension. 

3.1.1 Properties of the Alternative Bivariate 
Paradigm 

The differential expression 

{f^ + 2f^+f^)dxdy 

is readily seen to be rotation invariant, as is to be 
expected, given its physical meaning as the integrand 
for the representation (28) of the thin plate energy. It 
follows that every differential expression of the form 
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is also rotation invariant. This leads to the 



For different bivariate functions, however, it is pos- 
sible to establish non-uniqueness of a minimum. 
Consider a function 



Observation G: The Alternate Bivariate Functional 

e(f) is invariant under a rotation in the x, y-plane. 

L^t/=/(x, y) , g = g{x, j) G 5' be two interpolating 
spline functions. Invoking the general vector inequality 

||w|| + ||v|| > \\u + v||, we find 



^fl+^fl+fl+^8l+2gl+g. 



f +R 

J XX oxy 



+ 2 



( f +g t ( f + 

J xy oxv . J yy 



SO that 



2JQ 



'f^ + 



XX <JXX 



+2 



r^+; 



fn+< 



and in terms of the spline fiinctionals, 



e{f) + e{g) ^Jf+g 

2 2 



Since S is an affine space, the mean (/*+ g)l2 is again in 
S. Based on the above inequality we thus have 

Observation H: The alternate Bivariate functional 
(30) is a convex functional. In particular, as expressed 
in terms of the partial slopes, 

eif) = e(zxy zy^, ..., zx„, zyj, 

it is a convex function in terms of these partial slopes. 

3.1.2 A Related Non-Uniqueness Result 

We do not know at this point, whether the convexity 
is moreover strict. We expect that it is not. Strict con- 
vexity would imply uniqueness of minima and, again, 
we do not believe that this is the case. However, an 
example for such non-uniqueness is eluding us for our 
specific space of rHTC spline functions. 



s(x,y),-\ <x<+l,0<y 1 

representing a rectangular strip bent only in the direc- 
tion of the X-axis. In the direction of the j-axis, the 
function s(x, y) is constant (Fig. 15), in particular, for 
all J, 

s(-{,y)=-l,s(0,y) = 0,s(+Uy)=-l. 

The function s(x, y) thus interpolates the values -1,0, 
-1 along the contour lines in the direction of thej^-axis 
at the arguments x = -1 , 0, +1 . Given any Lavery spline 
f(x) considered in Sec. 2.3.3 as interpolating the exam- 
ple (21). Then it is readily seen, that any function 

minimizes both bivariate functionals: e(s), e(s). For 
functions corresponding to strips bent in only the strip 
direction, these bivariate functionals have thus multiple 
minima. This makes a strong case for non-uniqueness 
bivariate splines minimizing the functionals (29) and 
(30). However, we were not able to pinpoint non- 
uniqueness for the rHTC based spline functions consid- 
ered in this section. In particular, we were surprised to 
realize that the function s(x, y) cannot be generated 
using rHTC elements. As this fact points to an impor- 
tant limitation of using rHTC elements, we formulate 
the 




Fig. 15. Non-uniqueness of a minimuni for a bivariate function. 
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Observation J: The surface function of a rectan- 
gular strip bent only in the strip direction cannot be 
represented by means ofrHTC elements, 

3.2 Computing Bivariate Splines Based on an 
Alternate Paradigm 

The proposed alternate functional (53) is still expen- 
sive to minimize. In analogy to the procedure described 
in Sec. 2.4.1, we propose to approximate that function- 
al as follows. Let 



denote the area of the triangle 4 of the triangulation. 
Then the integral of each partial spline function or ele- 
ment y^ =fk{x, y) is individually approximated. 



[ 4fL+^f^+fL^y 



^ll,(JL+2f^+f^dxdy 



I (fL +24 +f^dxdy. 



This approximation suggests an iterative procedure 
based on minimizing weighted energy expressions: 
Starting with the classic spline f^\x,y), a sequence of 
interpolating spline functions is generated in hopes to 
converge towards the alternate Lavery integral (30), 

r\x,y),p\x,y),...,p\x,y),... 
with associated partial integrals 

^(//'^) = I if^' +2C + C)^xAy k = \, ..., m 
At each step / = 0, 1, 2, ..., weights 



w, 



,(0 



E{fn 



,k = l,..., m, 



are introduced. Given the function/^^"^\ the subsequent 
function/^^ is determined as the solution to the follow- 
ing minimization problem (compare (25)): 

(31) minX^r J,//r +2/^' +/^')dxdy 



For the exploratory implementation reported here, the 
above minimization of a quadratic function in the par- 
tial slopes 



ZX; , Zy: , ZX/ , ZV.- , ZX; , ZV; . 

is again formulated as solving a linear system of equa- 
tions in these variables, and Gauss-Seidel iteration pro- 
vided a workable option. Just as discussed in the uni- 
variate case, there are two iterative processes, one "on 
top of of the other. The "inner iteration" aims to solve 
the linear system of equations in order to solve the 
intermediate minimization problem (31), while the 
"outer" iteration sequentially creates spline functions 
expected to converge toward a limit that approximates 
the alternate Lavery spline defined by the minimization 
of the alternate Lavery functional (30). 

This approach raises the same issues as its parallel in 
the univariate case. The above iterative approximation 
based on consecutive minimization of weighted sums 
of integrals (3 1 ) cannot be expected to converge exact- 
ly to the alternate bivariate Lavery spline. The argu- 
ment relies on Observation F pertaining to the univari- 
ate case. Each iterate/^^(x, y) minimizes the energy of 
a physical structure consisting of thin plates of various 
thicknesses, namely the weights. Any minimizing 
shape of that structure, however, is known to have some 
vanishing second derivatives at its boundary. 
Univariate Lavery splines do not have this property, 
and we do not expect bivariate Lavery splines to have 
this property, either. 

Provisions need to be made should one of the inte- 
grals E{f^^~^^) vanish. This is a serious problem because 
very large or "infinite" weights will cause certain rHTC 
elements to remain stuck in linear form. In addition, 
limits of bivariate Lavery splines, in the original as well 
as the alternate formulation, may not be unique (see 
Sec. 3.1). These issues call for regularization proce- 
dures which are adequate for particular kinds of appli- 
cations. 

Indeed, spurious spikes were observed for the appli- 
cations demonstrated here. The problem was solved by 
interspersing an averaging step, where the partial 
slopes were replaced by averages of neighboring 
slopes. 

While our experiences with approximating alternate 
bivariate Lavery splines (Fig. 1 6) were encouraging, in 
particular, for very large data sets, much work needs to 
be done. 

The accuracy performance of the Gauss-Seidel itera- 
tion is a major issue. Simple examples show that initial 
steps are not even contracting. When does contraction 
start? Does it keep contracting after that? Much 
research has been aimed at these questions. These 
analyses of the convergence properties of the Gauss- 
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Iteration 1 




Iteration 3 



Iteration 5, Final 




fl fl 



Fig. 16. Example of 2-D Lavery splines. 



Seidel iteration, however, provide estimates of the 
number of iterations required to solve hnear equations 
to a specified accuracy which are much too large to be 
practical. To be practical, the number of "passes" 
through the entire system of linear equations on the 
order of 10 or 20, has to be sufficient. 

In our experience, relatively few steps of the "outer" 
iteration were sufficient to provide qualitatively satis- 
factory results. However, we do not have formal proof 
of convergence. 

Finally, and perhaps most importantly, a better 
understanding of the necessary regularization devices 
needs to be developed. 



4. Conclusions 

In this paper we have investigated some properties of 
non-oscillatory splines introduced by John Lavery 
[11,12]. These splines, called in this paper Lavery 
splines, minimize what we have termed the Lavery 
integral (4), We have seen that the minimizing spline 
for (4) does indeed model sharp edges and jumps in 
data without introducing the "Gibbs phenomenon" at 
the comers. We have shown that the Lavery integral 
and the associated Lavery splines satisfy a number of 
properties. First, we have shown that the Lavery inte- 
gral of a cubic spline function is the absolute value of 
an alternating sum of inflection slopes and of the end 



slopes. Next, we showed that any positive linear com- 
bination of Lavery splines for the same interpolation 
problem is again a Lavery spline for the same problem. 
Furthermore, two Lavery splines for the same opti- 
mization problem share essentially the same inflection 
points and that the free ends of Lavery splines are either 
linear fimctions or they contain an inflection point in 
the interval of definition. 

Two algorithms for estimating Lavery splines have 
also been considered. The first algorithm introduced by 
John Lavery [11,12] reduces to solving a least absolute 
value minimization problem for which he used an inte- 
rior point method for linear programming to obtain the 
minimum spline coefficients. The absolute value mini- 
mization was based on a "discretization" which leads to 
a very large number of variables as sufficiently smaller 
discretizations were considered. The extension of this 
method to bivariate Lavery splines also led to compu- 
tationally intensive compute times even for moderate 
data sizes. The authors have introduced a modified 
approach based on an iterated weighted least squares 
algorithm. Although the minimizing spline this algo- 
rithm produces is not a Lavery spline it is an approxi- 
mation that also produces sharp edges without the 
"Gibbs phenomenon," 

In the bivariate case the Lavery integral minimizes 
the integral of the sum of the absolute values of the sec- 
ond partial derivatives of the spline. Analogous to the 
univariate case we introduce an alternative variational 
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principle based on the integral of the square root of the 
sum of squares of the second partial derivatives of the 
bivariate spline. Whereas the Lavery integral in the 
bivariate case is not rotationally invariant, the alterna- 
tive principle is. The alternative bivariate functional is 
shown to be a convex functional and, when expressed 
in terms of the partial slopes at the vertices of an under- 
lying element, it is a convex function of these slopes. 
We have also shown that a rectangular strip bent only 
in the strip direction cannot be represented by means of 
a triangular rHTC element. 

Analogous to the alternative weighted least squares 
algorithm introduced to compute the approximate uni- 
variate splines we introduce an alternative weighted 
least squares algorithm in the bivariate case. Although 
we have not introduced a convergence analysis of the 
alternative algorithms in this paper, computational 
experience has shown that the iterative, Gauss-Seidel 
based, algorithm used provides qualitatively satisfacto- 
ry results in only a few iterations of the main loop. 
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